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i  . 

! A  splitting-up  method  and  an  implicit  finite  difference 
method  are  presented  to  solve  time-dependent,  one-dimensional, 
laminar,  premixed  flame  problems.  An  example  for  studying  the 
development  of  an  ozone  decomposition  flame  is  calculated.  A 
movable  boundary  technique  is  adopted,  therefore  the  grid 
points  can  be  significantly  reduced.  Special  care  is  taken  to 
maintain  the  accuracy  of  the  solution.  The  results  are  checked 
in  many  ways.  All  checks  show  that  the  present  method  is 
satisfactory.  . 
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SIGNIFICANCE  AND  EXPLANATION 


More  and  more  attention  is  paid  to  combustion  problems 
not  only  by  engineers  but  also  by  mathematicians,  because  a 
number  of  interesting  and  difficult  problems  occur.  For 
example,  a  premixed  flame  problem  will  reduce  to  a  typical 
reaction-dif fusion  equation. 

It  is  well  known  that  in  a  premixed  combustible  fluid 
mixture  a  steady  flame  will  be  developed  when  it  is  ignited. 
This  fact  has  been  proved  theoretically  for  a  simple  chemical 
reaction  model.  There  have  been  many  works  on  studying  the 
conf iguration  of  a  laminar,  premixed  flame.  These  works  have 
essentially  followed  two  different  approaches.  One  is  the 
steady-state  approach  in  which  the  problem  to  be  solved  is 
reduced  to  a  two-point  boundary  value  problem  of  ordinary 
differential  equations.  The  other  is  the  time-dependent 
approach  in  which  the  full  time-dependent  equations  with  proper 
boundary  and  initial  conditions  are  treated.  In  this  paper 
our  attention  is  focused  on  the  latter  approach. 


Since  the  combustion  processes  are  highly  exothermic, 
chemically  reaction  processes  and  there  exists  a  number  of 
vastly  differing  time  scales,  numerical  solutions  will  suffer 
from  a  number  of  difficulties.  One  is  stiffness.  With  this 
in  mind  a  splitting-up  method  is  presented,  which  splits  the 
chemical  kinetic  terms  from  the  fluid  mechanical  terms.  We 
think  it  might  ameliorate  some  of  the  difficulties.  At  the 
same  time,  an  implicit  finite  difference  method  is  presented 
in  order  to  compare  and  identify  the  validity  of  the  methods. 

Generally  speaking,  the  region  of  calculation  must  be 
taken  large  enough,  therefore  a  large  mount  of  grid  points 
must  be  taken.  Obviously  it  is  costly.  In  order  to  reduce 
grid  points  a  movable  boundary  technique  is  adopted. 

Special  care  of  error  control  is  taken  to  maintain  the 
accuracy  of  the  solutions. 


For  comparison  of  the  results  obtained  by  the  present 
methods  with  published  results  an  example  for  an  ozone 
decomposition  flame  is  calculated.  The  results  are  also 
checked  in  many  ways.  The  comparison  and  check  .siiaw.-that-t-h' 
present  methods  are  satisfactory.  P  \j  | 
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The  responsibility  for  the  wording  and  views  expressed  in  this 
descriptive  summary  lies  with  MRC,  and  not  with  the  author  of 
this  report. 


NOMENCLATURE 


space  coordinate 
time  coordinate 
pressure 

density  of  fluid  mixture 
temperature 

mass  fraction  of  the  i-th  species 
velocity  of  fluid  mixture 
universal  gas  constant 

specific  heat  capacity  at  constant  pressure  of 
the  fluid  mixture 

specific  heat  capacity  at  constant  pressure  of 
the  i-th  species 

rate  of  production  of  i-th  species 
molecular  weight  of  i-th  species 

binary  diffusion  coefficient  for  the  i-th  species 
thermal  conductivity 
enthalpy  of  i-th  species 

standard  enthalpy  of  formation  of  i-th  species 


NUMERICAL  METHODS  FOR  SOLVING  A  PREMIXED  LAMINAR  FLAME 

Xi-Chang  Zhong 

1.  Introduction 

More  and  more  attention  is  paid  to  combustion  problems 
not  only  by  engineers  but  also  by  mathematicians,  because 
a  number  of  interesting  and  difficult  problems  occur.  For 
example,  a  premixed  flame  problem  will  reduce  to  a  typical 
reaction-diffusion  equation. 

It  is  well  known  that  in  a  premixed  combustible  fluid 
mixture  a  steady  flame  will  be  developed  when  it  is  ignited. 
This  fact  has  been  proved  theoretically  for  a  simple  chemical 
reaction  model.  There  have  been  many  works  on  studying  the 
configuration  of  a  laminar,  premixed  flame.  These  works 
have  essentially  followed  two  different  approaches.  One  is 
the  steady-state  approach  in  which  the  problem  to  be  solved 
is  reduced  to  a  two-point  boundary  value  problem  of  ordinary 
differential  equations.  The  other  is  the  time-dependent 
approach  in  which  the  full  time-dependent  equations  with 
proper  boundary  and  initial  conditions  are  treated.  In  this 
paper  our  attention  is  focused  on  the  latter  approach. 

Since  the  combustion  processes  are  highly  exothermic, 
chemically  reaction  processes  and  there  exists  a  number 
of  vastly  differing  time  scales,  numerical  solutions 
will  suffer  from  a  number  of  difficulties.  One  is 
stiffness.  With  this  in  mind  a  splitting-up  method  is 
presented,  which  splits  the  chemical  kinetic  terms  from  the 

Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29 
80-C-004 1 . 


fluid  mechanical  terms.  We  think  it  might  ameliorate  some 
of  the  difficulties.  At  the  same  time,  an  implicit  finite 
difference  method  is  presented  in  order  to  compare  and  identify 
the  validity  of  the  methods. 

Generally  speaking,  the  region  of  calculation  must  be  taken 
large  enough,  therefore  a  large  mount  of  grid  points  must 
be  taken.  Obviously  it  is  costly.  In  order  to  reduce  grid 
points  a  movable  boundary  technique  is  adopted. 

Special  care  of  error  control  is  taken  to  maintain  the 
accuracy  of  the  solutions. 

For  comparison  of  the  results  obtained  by  the  present 
methods  with  published  results  an  example  for  an  ozone 
decomposition  flame  is  calculated.  The  results  are  also 
checked  in  many  ways.  The  comparison  and  check  show  that 
the  present  methods  are  satisfactory. 


2 


Formulation  of  the  Problem 


2.1  Governing  Equations 

We  consider  one-dimensional  flow  and  neglect  the  effects 
of  radiative  heat  transfer  and  thermal  and  pressure  diffusion. 
The  governing  equations  are  as  follows. 

Continuity 


3P  .  3(Pu)  _  „ 

St  +  Sx -  "  0  ' 


(2.1) 


Conservation  of  momentum 


p  +  pu  =  _3R+J_  r{u  +  iK) 

M  Jt  u  9x  9x  9x  liy  3  '  9xJ  ' 


(2.2) 


Conservation  of  species 


3Yi  3Yi 

P  StT  +  Pu  asr 


3y< 


"  &  <p  Di  S5T>  +  “i  • 


(2.3) 


Conservation  of  energy 


PC  +  PC  u  |^  =  +  u  +  (y+  £  k)  (|H.)2 

p  3t  p  9x  3t  9x  3  9x 


+  -  <>  -  £  Vi 


55 

N 


9yi  9T 


-  Ji  PDtCPi  ST  55  - 


(2.4) 


and  equation  of  state 


P  = 


mH  RT 


(2.5) 


The  variables  apearing  in  these  equations  have  their  usual 


meaning  as  listed  in  nomenclature.  In  order  to  simplify 
these  equations  further,  the  effect  of  viscosity  is  assumed 
to  be  negligible  and  fluid  velocity  is  assumed  to  be  small 
compared  to  the  speed  of  sound.  From  the  latter  we  can 
integrate  equation  (2.2)  and  obtain  the  following  condition. 

p  =  P0  =  const.  (2. 6) 

Then  the  equations  (2.1)  -  (2.5)  are  simplified  into 
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where  mn  =  Pu  .  Since 
u  '  x=0 


3  _  3^  3  _3_ 

Tx  3x  Tip-  “  P  3<j> 

3  34»  3  .  3  .  _  .  .3^3 

3t  "  3t  +  3t  =  Pu  +  m0)  W  +  %t  ' 


the  remaining  equations  become 

3yi  3yi  a  2  3yi  “i 

5t  ~  +  m0  W  =  5F  (P  Di  ^T")  +  T 


(2.12) 


3t  a.  m  8T  -  1  9  3t\  ?  M  K 

Tt  +  m0  5UT  ~  c~  W  (PX  W}  ~  A  Wihi 

p  i=l 


N  p.  ~ 

-  *  c^pS 

i=l  cp  1 


3yi  3T 

i  "51J5-  DJT 


(2.13) 


By  introducing  the  following  nondimensional  variables: 
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the  equations  can  be  nondimensionalized  as 


9yi 


+  m0  =  <L<J  1  (P2Di  (2.14) 


3yi. 


u) . 


3T  .  3T 
Tt  +  m0  W 


<r  *  '°x  H>  -  re-  ?  Vi 

P  P  1=1 


C 

N  V 


3y< 


-  r1  y  ^  p2d  1  8T 
(L«)  ,£,  ~  p  Di  SF"  ^  ' 

1-1  p 

x 


where 


L  = 


00  P~DC 

00  "  P„ 


(2.15) 


is  a  characteristic  Lewis  number,  and  P  ,  T  ,  D  ,  X  ,  c 

'  OO  OO  OO  '  CO  n 

^  00 

and  are  some  characteristic  values  of  density,  tempera¬ 
ture,  mass  diffusivity,  heat  conductivity,  specific  heat 
capacity,  and  molecular  weight,  respectively.  For  simplicity, 
the  superscript  star  is  omitted  in  the  above  equations. 


2.2  Initial  and  boundary  conditions 

The  equations  (2. 14) - (2 .15)  are  a  parabolic  system; 
thus  boundary  and  initial  conditions  are  required.  In 
this  paper  we  consider  the  propagation  of  a  premixed  flame. 
So  we  can  specify  the  boundary  conditions  in  the  following 
way.  At  the  burned  boundary  the  burned  values  are  taken 
and  at  the  unburned  boundary  the  unburned  values  are  taken, 
i.e. 
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T  =  Tfa  / 

at  burned  boundary 

yi  =  *i  b  ' 

T  =  Tu  , 

at  unburned  boundary 

yi  =  y.  • 

u 

The  unburned  and  burned  values  must  satisfy  some  condi¬ 
tions,  which  depend  on  what  model  is  assumed.  For  example, 
in  this  paper  we  assume  that  the  flame  is  an  adiabatic 
flame  and  the  burned  values  satisfy  chemically  equilibrium 
equations . 

The  initial  conditions  are  specified  as  a  step  function. 


f  T 

¥ 

> 

* 

T(*,0)  = 

1  u 

i 

for 

c 

i 

l  Tb 

¥ 

< 

*c 

yi(^*0)  =  | 

r  yiu 

for 

¥ 

> 

*c 

1  yib 

>P 

< 

*c 

where  <l<  is  a  given  value, 
c 


(2.16) 
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2.3  Chemical  Kinetic  Formulation 

The  elementary  chemical  reactions  are  expressed  as 


N 

I 

i=l 


m,  x 


Xi 


N 

l  v"  . 
i=l  m'1 


Xi 


m  =  1 ,2,  . . . ,M.  (2.17) 


The  rate  of  production  of  the  i-th  species  appearing  in 
equations  (2.12),  (2.13)  is  given  by  the  law  of  mass  action. 


u . 
i 


M 

l 


m=l 


(v"  < 
m,  i 


(2.18) 


where  v'  ,  and  v"  .  are  the  stochastic  coefficients  of  the 
m,  x  m,  x 

species  i,  i  =  1,...,N,  appearing  in  a  reactant  and  product 
in  the  reversible  reaction  m,  m  =  1,2,...,M.  The  c^  are 
the  moles  per  unit  volume  of  species  i  and  related  to  the  mass 
fraction  by 


c  =  JL  y 

i  Mi  yi 


(2.19) 


The  specific  rate  constants  for  forward  and  backward  mode 
of  reaction  m  are  usually  given  by  the  following  expression. 

*  f  Si  “Em/RT 
k*  =  Bf  T  m  e 
m  m 


(2.20) 


B 


m 


f  f 

,  s_  ,  E_  =  const. 


m 


m 


b  b 

has  a  similar  expression.  The  constants  e£  and  E„  are 
m  mm 

the  activation  energy  of  the  forward  and  backward  mode  of 


lb  iiiw  'i 
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reaction  m,  respectively.  In  terms  of  nondimens ional 
variables  hk  in  the  equation  (2.14)  is  expressed  as 


w.  M 


— 5-  =  I  M.  <v"  .  -  V'  . )  [k*  P 
«  •  .  i  m  / 1  m  ,  x  m 

P  i=l 


a  -1  , 

*  *  m  ,N  y.  v*  . 
1 — r  t _ *  \ 


TJ 

1=1  i 


*  *em“1  N  y.  v"  . 

p  rr  (^-)  m,:L 


'  *  M. 

1=1  1 


(2.21) 


where  the  superscript  star  refers  to  the  nondimens ional 
variable  and 


3  =  l  v "  . 

m  i^1  m , l 


a  =  l  v'  .  , 
m  m,  l 


*  p  a  -1 

kf  =  kf  m 

m  m  <*> 


*  P  3  -1 

kb  =  kb  m 

m  m  °°  M_ 


Correspondingly,  in  the  equation  (2.14) 


1  N  *  *  M  N  hj 

~  y  h,  w .  =  i  i  ..  —  + 


(2.22) 


—  y  h  w.  =  y  y  < - - - +  c  (t  -t")my(v"  .-v*  .)) 

*  .4.  i  i  '  M . C  T  p.  0'"i'  m, l  m,i" 


P  i=l 


x=l  i  P™  "  x 

*  *am_1  N  y.  v'  .  K*  .Pm-1  N  y.  v"  . 

*  K  p  TT  (m1)  m'1-  km  p  TT*^)  m,1i 

i=l  i  i=l  1 


(2.23) 


3. 


Numerical  Method 


3.1  Coordinate  Transformation 

Strictly  speaking,  the  infinite  boundary  must  be 
considered,  but  it  is  impossible  and  unnecessary  in  practice. 
In  general  one  can  take  a  value,  4>  large  enough  to  allow  the 
full  development  of  the  flame  before  any  effect  could  be 
felt  there.  However  it  is  costly  to  do  so.  In  order  to 
reduce  grid  points,  it  is  effective  to  adopt  the  movable 
boundary  and  to  introduce  the  following  coordinate  trans¬ 
formation: 

♦  -  ¥b(t)  ¥  -  Yb(t) 

?  =  ^u(t)  -  4>b(t)  =  e  ( t) 

t  =  t 

where  ^  ,  <1*  represent  the  burned  and  unburned  values, 

,  respectively.  Since 

3  _  3  3£  .  3 

TE  ~  K  "5t  +  TE 

3  _  1  3 

e  "5T 

the  equations  (2.14) 
m0+b  3yi 

W  +  ~  W 

3T  ,  m0+b  3T 
3t  e  3T" 


-  (2.15)  become 


3  roV  yii  +  i 

W  (p  D1  9X-5  +  T 


(3.1) 


^  n  <pA  $  “  ps-  £  ^ihi 

Cpe  p  i-l 

1  N  Cp.  3y.  3T 

'  r?  Ji  c7" Di  (3'2’ 
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where  b  =  -  -  $.)  +  and  the  dot  represents  the  deriva- 

u  b  D 

tive  with  respect  to  t. 

In  order  to  solve  the  equations  (3.1) -(3. 2),  two  numeri¬ 
cal  methods  are  adopted  here.  One  is  the  splitting-up  method 
and  the  other  is  the  implicit  finite  difference  method.  In 
what  follows,  we  describe  these  methods  respectively. 

3.2  Splitting-up  Method 

The  system  of  equations  (3.1) -(3. 2)  is  usually  stiff 
and  the  solution  frequently  has  rapid  transients.  In 
other  words,  there  are  a  number  of  time  scales  and  these 
scales  vary  widely.  For  example,  in  the  vicinity  of  the 
flame,  the  chemical  reactions  are  quite  rapid,  compared 
with  the  fluid  mechanical  scale  and  among  the  chemical 
reactions  some  may  be  quite  faster  than  others.  The 
splitting-up  method  could  allow  one  to  use  different  step 
sizes  according  to  different  time  scales,  e.g.  one  large 
step  size  can  be  used  to  advance  the  slowly  varying  terms 
while  several  smaller  step  sizes  can  be  taken  on  the  faster 
terms.  For  simplicity  we  rewrite  the  equations  (3.1) -(3. 2) 
as 

|~=a^-+d^-(n||)  +  Q  (3.3) 

We  split  equation  (3.3)  into  two  parts  which  group 
the  fluid  mechanics  and  echmical  kinetics  terms: 
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I 


Obviously,  the  equation  (3.6)  is  a  tridiagonal  system. 

So  the  solutions  for  any  f  and  all  j  can  be  accomplished 
with  a  tridiagonal  matrix  algorithm. 

Having  obtained  the  solution  for  the  first  part  of 
the  split,  we  may  turn  our  attention  to  the  equation  (3.5). 
The  solutions  f can  be  thought  of  as  predicated 
intermediate  values  of  the  solution,  and  are  used  as  the 
initial  conditions  of  the  equation  (3.5).  Since  g  contains 
no  spatial  derivatives  equation  (3.5)  is  a  system  of 
ordinary  differential  equations  at  each  grid  j.  Usually 
the  equation  is  stiff,  so  some  effective  stiff  methods 
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such  as  the  Hindmarsh-Gear  package  can  be  used. 

In  order  to  improve  the  accuracy  and  efficiency  of 
the  calculation,  a  symmetric  split  operator  is  used. 

Let  LD  be  the  operator  providing  a  predicated  value 
of  f and  be  the  operator  that  obtains  a  correc¬ 
tion  of  that  first  step.  The  following  operator  is  used 
to  advance  the  solution  from  t  to  t  +  2At: 


f  i+2 

j 


-  LDLCLCLDfj 


3.3  Implicit  Finite  Difference  Method 

In  order  to  compare  and  identify  the  validity  of 
the  methods,  an  implicit  finite  difference  method  is  also 
presented.  For  equation  (3.3)  the  implicit  finite 
difference  scheme  is  as  follows. 


a .  i+1  i+1  d .  i+1 

=  id  (fj+l  "  fj  -l)  +  ^2  [nj+l/2fj+l 

i+1  i+1 

-  (nj+l/2+nj-1/2)  fj  +nj-l/2fj-l]  + 


(3 


where  the  index  has  the  same  meaning  as  in  (3.6). 
Similarly,  the  tridiagonal  matrix  algorithm  is  also  used. 


3.4  Error  Control 


In  calculating  the  nonlinear  terms  a,  d,  n,  g,  of 
equations  (3.6)  and  (3.8)  are  linearized  and  evaluated 
explicitly.  In  order  to  control  error  the  following 
approach  is  taken.  We  take  two  steps  of  At/2  and  one 
step  of  At;  then  compare  the  solutions  at  t+At.  Let 


e 


max 

j 


pji+1  (1) 
max 


_  ji+1  (2) 

|ij+1  (2,l  _ 


(3.9) 


where  superscripts  (1),  (2)  refer  to  the  number  of  steps. 

Whether  the  step  size  is  increased  or  decreased  depends  on 
whether  e  is  less  than  or  greater  than  E.  E  is  an  error 
tolerance  and  is  specified  beforehand.  If  the  convergence 
has  been  reached  (i.e.  e  <  E) ,  ^  is  taken  as  the 

solution  at  t  +  At. 

3.5  Movable  Boundary  Technique 

As  previously  pointed  out,  if  the  boundary  is  fixed, 
a  large  value  of  <P  must  be  taken.  Thus  it  is  computationally 
costly.  It  is  naturally  desirable  that  the  calculated  region 
be  confined  to  such  that  it  always  contains  the 
flame  and  it  is  as  small  as  possible.  It  is  well 
known  that  in  a  premixed  combustible  fluid  mixture  a  steady 
flame  will  be  established  after  ignition.  In  other  words, 
the  flame  will  propagate  through  the  combustible  fluid  with 
a  constant  velocity.  In  view  of  this  fact,  the  movable 
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boundaries  could  be  used.  Moreover,  the  moving  velocity 

could  be  taken  to  be  the  same  as  the  velocity  of  the  flame. 

•  • 

Let  ^u,  'I'k  ,  Su  represent  the  moving  velocities  of  the 
unburned  boundary,  the  burned  boundary,  and  the  flame, 
respectively.  Then  the  boundary  values  of  ^  are  obtained 
by  integrating 


d*u 

dt 


d*b 

dt 


(3.10) 


The  flame  velocity  based  on  the  density  of  the  unburned  mixture 
can  be  evaluated  by  means  of  the  mass  conservation  principle. 
That  is,  the  net  mass  production  rate  of  species  i,  inside 
a  control  volume  (which  is  large  compared  to  the  thickness 
of  the  flame) ,  must  be  eoual  to  the  mass  rate  of  outflow 
of  species  i.  The  resulting  equation  in  terms  of  is 

b  ». 

-p1  /  pu(yib-nu>  (3-n) 

u 


In  principle,  any  of  the  chemical  species  can  be  used 
to  compute  the  flame  velocity,  because  all  these  computed 
flame  velocities  should  be  identical.  But  in  practice, 
owing  to  the  numerical  error  there  are  some  differences 
between  these  velocities.  In  our  calculation  one  of  them 
is  taken  and  the  other  is  used  to  check  the  accuracy  of  the 
calculation. 
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In  order  that  the  calculated  region  can  be  automatically 
confined  to  where  significant  changes  occur,  the  grid 
interval  is  allowed  to  expand  or  contract.  The  way  to  do  so 
is  as  follows.  At  the  beginning,  we  use 


•  • 

Ip  =  Ip  , 
TU  U1 

y.  -  y 

J1U  1 

XT  w  ' 

L^iu 

l  =  l  . 

vui 

2  -  {( 

' 1 

~  t 

>b  ■». 

ib 


^ibb 


i  )  x  d 


lu 


>1 


(3.12) 


where  <Jiui  =  -J^  dip  /  (yib  -  yiu>  ,  and  6, a  are  some 
constants  (e.g.  6  =  100,  a  =  0.1).  Subscripts  uu,  bb 
represent  some  grid  points  near  the  unburned  and  burned 
boundaries,  respectively.  After  the  calculation  proceeds 
some  time  we  turn  to  use  y. 

r  b 

L  a. 


K-K-  •pusu 


(3.13) 


yib-y 


iu 


The  reason  for  doing  so  is  that  we  found  it  was  not  easy 
for  the  flame  velocity  to  be  stable  if  the  formulas  (3.12) 
are  used  only. 


-16- 


4 


Example  of  Calculation 


The  methods  described  above  are  applied  to  calculate 
the  structure  of  an  ozone  decomposition  flame.  In  order 
to  compare  the  results  with  the  published  results,  the 
following  nonessential  approximations  were  made: 

2  2 

0^=02=  ...  =  =  D  and  P  D  =  const.  =  , 

pX  =  const.  =  P  X 

CO  oo 


The  reaction  mechanism  for  the  ozone?  ^r.o;v;©s ition  consists 
of  seven  reversible  reactions: 

02  +  0  +  X  (4.1) 


202  (4.2) 

20  +  X  (4.3) 

Where  X  represents  any  of  the  three  species  0,  02#  O^- 
Below  we  illustrate  how  to  determine  the  boundary 
conditions.  At  unburned  boundary  a  combustible  mixture 
of  75%  O2  and  25%  0^  (by  volume)  at  a  temperature  of  300°K 
are  assumed.  That  is. 


0,  +  X 
3  <vT 


0+0, 


k4 

f 


0-  +  X 
2  XT 


*lu  =  0 


Tu  =  300 °K 


(4.4) 


The  burned  values  can  be  determined  by  the  unburned  values. 
Between  these  values  some  conditions  must  be  satisfied. 

In  this  case,  they  are 

Conservation  of  total  enthalphy 

C,,b-Cp,.‘!lbh!1'=>3.lll!l  <4-5 

Chemically  equilibrium  equation 


The  equations  (4.5),  (4.6)  are  a  system  of  nonlinear  equations.  The 
Newton  iteration  method  is  used.  The  solutions  are 

y..  =  0.1259  x  io~7  ,  y  =  l  -  0. 1259x10-7 , 
lb  2b  (4.7) 

y3b  "  0  '  Tb  =  1246. 9»K. 

N 

The  rate  of  production  w.  and  the  term  (  l  h.<u.)/p 

1  i=l  1 

are  given  by  the  equations  (2.21)  and  (2.22).  The 
thermodynamics  and  kinetics  data  used  here  are  given 
in  Table  I. 


The  initial  conditions  were  taken  to  be 


TABLE  I 


DATA  FOR  OZONE  DECOMPOSITION 


Symbol 

Value 

Symbol 

Value 

Ef'EI-Ef 

24140  -isffl 

Bb 

B4 

1.88  x 

106 

E4 

0000  -isfl 

B5*  B6 ,B7 

2.47  x 

102 

E5'E6'E7 

117350  mo?e 

u(D 

ho 

58675 

cal 

mole 

0 

u  ( 2) 

ho 

0 

E4 

95“10  mole 

K  ( 3) 
h0 

34535 

cal 

mole 

E5'E6'E7 

0 

p» 

1.201X10"3  -3-j 
cnr 

f  f  f 
^1 ' s3 

5/2 

T 

OO 

300. 

O 

K 

f 

5/2 

0.2524 

cal 

S4 

CPoo 

g-°K 

sf  sf  sf 

S5'S6'S7 

5/2 

9.112X10'5 

cm- sec 

b  b  b 

7/2 

n 

A 

OO 

sl's2's3 

OO 

Pc  Le 

cd  n  oo 

“oo 

c,b 

5/2 

Moo 

16 

_  9,  . 

S4 

mole 

s5's6'S7 

7/2 

4 . 203x10~3  cm 

6.76  x  106 

5.878X10-5  sec. 

Bf 

B4 

4.58  x  106 

*»/*«, 

71.51 

cm 

sec 

Bf  Bf  Bf 
B5'  6'  7 

5.71 x  106 

po 

0.821 

atm. 

B1-B2'B3 

1.18  x  102 

M 

16 

9 

Mi 

mole 

*2 

32 

9 

mole 

M3 

48 

9 

mole 
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i 

i 


iksthsLx^-vi  ~ 


(4.8) 


y^io 

y2(*) 


y3  (^) 


! 


0 

0.1259  x  10"7 

2 

I 

1-  0.1259*10~7 

1 

T 

0 


*  >  *c 

*  I  *c 

*  >  *c 

*  I  *c 

♦  >  *c 
*  I  *c 


T(IJ>) 


300  °K 
.  1246. 9°K 


*  >  *c 

*  i  *c 


(4.9) 


(4.10) 


(4.11) 


5.  Results  and  Discussion 

The  solution  of  the  equations  (3.1) -(3. 2)  with 
boundary  conditions  (4. 4) -(4. 7)  and  initial  conditions 
(4 . 8) - (4 . 11)  have  been  obtained  for  L_  =  1.  In  the 
following  we  describe  some  results  and  discuss  them. 

In  order  to  compare  and  identify  the  validity  of  the 
methods,  the  splitting-up  and  implicit  finite  difference 
method  are  used.  In  the  former  the  integration  of  ODE  (3.5) 
uses  the  DGEAR  routine  of  the  IMSL  package,  which  is  an 
adaptation  of  the  package  designed  by  Hindmarsh-Gear .  The 
steady  profiles  of  temperature  and  concentrations  for  a  fully 
developed  ozone  flame  are  shown  in  Figures  1-2.  Here  29 
grid  points  are  used  and  =  12.97.  From  the  figures 

it  is  seen  that  the  agreement  of  the  results  obtained  by  two 
methods  is  good.  The  flame  velocities  based  on  O2  are 
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aSSitea^iBiaiga 


implicit 


f 


Su  =  51.5  cm/sec  for  the  splitting  method, 

Su  =  49.7  cm/sec  for  the  implicit  method. 

Margolis  [ 1 ]  used  the  method  of  lines  involving  collocation 
with  B-splines  and  obtained  Su  =  49.7  cm/sec,  while 
Bledjian  [2]  used  the  method  of  lines  with  low  order  accuracy 
and  got  Su  =  54  cm/sec. 

We  have  analyzed  the  eigenvalues  of  the  Jacobian 
derivative  matrix  for  the  chemical  reaction  equation 

dyi  .  “i 
dt  P  * 

For  some  typical  cases,  e.g.  the  fully  burned  mixture,  we 
obtained  the  ratio  of  maximum  and  minimum  eigenvalue 
is  about  3*10  .  That  means  the  stiffness  is  mild  in  the 
present  example.  So  it  is  not  surprising  that  two  methods 
can  obtain  satisfactory  results.  However,  we  think  that 
if  the  stiffness  becomes  serious,  the  splitting-up  method 
might  be  more  effective. 

In  addition,  the  flame  velocity  is  senstive  to  the 
accuracy  of  variables,  so  the  appropriate  control  error 

must  be  chosen.  In  the  present  case,  the  control  error 

-4  -5  . 

E  is  taken  to  be  10  in  equation  (3.9)  and  10  in 

integration  of  the  ODE. 

In  order  to  establish  the  validity  of  the  movable 
boundary  technique,  the  fixed  boundary  is  also  considered. 

We  take  ||»  -  =  50  and  99  grid  points.  The  time  development 

of  the  right  propagating  flame  and  the  profiles  of 


temperature  and  concentrations  are  shown  in  Figures  3-5. 

The  flame  velocity  is  in  good  agreement  with  that  obtained 
by  movable  boundaries.  Their  deviation  is  about  0.7%. 

From  this  it  is  seen  that  the  movable  boundary  technique 
can  significantly  reduce  grid  points;  therefore  saving 
machine  time. 

As  for  the  grid  points,  some  tests  are  made. 

When  grid  points  are  reduced  to  19,  the  results  still  remain 
good.  The  flame  velocity  is  almost  the  same  as  that  obtained 
by  29  grid  points  (the  difference  between  them  is  only  in 
the  fourth  decimal) .  However,  when  grid  points  are  reduced 
to  14,  the  results  get  worse.  The  results  with  19  grid  points 
are  shown  in  Figures  6-7. 

The  results  are  checked  in  many  ways.  It  is  well  known 
that  the  concentrations  must  satisfy 

N 

l  Yi  *  1  • 
i=l  1 

In  the  present  calculation  this  condition  is  satisfied  very 
well.  Its  deviation  from  1  is  only  in  the  fifth  decimal. 

As  mentioned  above,  a  comparison  between  the  flame 
velocities  based  on  different  species  can  be  taken  as  a  check 
on  the  results.  In  our  calculation  the  difference  between  the 
flame  velocities  based  on  and  0^  is  less  than  0.2%. 

Another  different  case  is  calculated,  too,  in  which 
the  density  P  is  considered  as  a  constant.  It  has  been 
proved  that  when  P  is  constant  there  also  exists  a  steady 
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flame  [5].  Our  calculation  shows  that  it  is  indeed  true. 

The  steady  profile  of  temperature  and  concentrations  for 
an  ozone  flame  with  constant  density  is  shown  in  Figures  8-9. 

Finally,  a  remark  has  to  be  made.  When  using  the 
splitting-up  method,  one  has  to  decide  whether  the  solution 
is  considered  after  an  sequence  or  after  an 

LcLdLdLc  sequence.  Generally  speaking,  both  yield  valid 
solutions,  but  there  are  always  some  differences  between  them. 
In  order  that  the  solutions  after  an  or  an  step 
are  within  a  certain  tolerance  of  each  other,  an  additional 
constraint  on  the  step  size  is  required.  Fortunately  when 
using  the  control  error  mentioned  above,  the  solution  after 
an  Lc  is  in  good  agreement  with  that  after  an  1^. 

The  difference  between  two  flame  velocities  is  less  than  0.2%. 

7.  Conclusion 

A  splitting-up  method  and  an  implicit  finite  difference 
method  have  been  presented.  An  example  for  studying  the 
development  of  an  ozone  decomposition  flame  is  calculated. 

The  calculation  shows  that  two  methods  can  obtain  good  results. 
But  the  splitting-up  method  should  be  more  effective  when  the 
stiffness  becomes  serious.  A  movable  boundary  technique 
is  adopted,  therefore  the  grid  points  can  be  significantly 
reduced.  The  results  are  checked  in  many  ways.  All  checks 
show  that  the  methods  are  satisfactory.  The  present  method 
can  be  applied  to  compute  more  general,  time-dependent. 
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one-dimensional,  laminar  flame  problems 
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